# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)library(gridExtra)library(png)# tablelibrary(gt)library(gtable)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar angle and timezonelibrary(suncalc)library(lutz)# data wranglinglibrary(tidyr)library(dplyr)library(data.table)library(purrr)library(forcats)library(lubridate)
2 Setting up custom function
2.1windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet_wrap = F,facet_grid = F,hour_sunset =20,hour_sunrise =7) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, hour_sunrise, hour_sunset), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin labelslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet_wrap =="origin") { data2$origin <-unique(data_to_plot$origin) } elseif (facet_wrap =="trip") { data2$trip <-unique(data_to_plot$trip) }# add trip and origin if anyif (facet_grid) { data2$origin <-unique(data_to_plot$origin) data2$trip <-unique(data_to_plot$trip) }# set up max value maxvalue <-20# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="black",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet_wrap =="origin") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) } elseif (facet_wrap =="trip") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(trip ~ .,strip.position ="right" ) }# if facetif (facet_grid) { plt.dirrose_2 <- plt.dirrose_2 +facet_grid2(trip ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
2.2matlab_to_posix
Code
# function to convert matlab date into R posixmatlab_to_posix <-function(x, timez ="UTC") { days <- x -719529# 719529 = days from 1-1-0000 to 1-1-1970 secs <- days *86400# 86400 seconds in a dayreturn(as.POSIXct(strftime(as.POSIXct(secs,origin ="1970-1-1",tz ="UTC" ),format ="%Y-%m-%d %H:%M:%S",tz ="UTC",usetz =FALSE ),tz = timez ))}
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")# remove id == `2016014` and id = `2018002` (no benthic dives)data_dive <- data_dive %>%filter(!id %in%c("2016014", "2018002"))
We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.
Code
compute_table_1_a <-function(x) {# if the sum of the column is below 100if (sum(x) <100) {# it's definitely a proportion# so we weighted average these y <-weighted.mean(x, c(1068746, 2207361)) y <-paste0(round(y *100, 1), "%")return(y)# otherwise we just sum the column } else { y <-prettyNum(sum(x), big.mark =",")return(y) }}
# custom function to summarize each column differentlycompute_table_1_b <-function(x) {# if the column is numericif (all(is.numeric(x))) {# and the sum is below 100if (sum(x) <100) {# it is time# so we "circular" weighted average them y <- psych::circadian.mean(c(rep(x[1], 237),rep(x[2], 187) ),na.rm = T ) y <-round(y, 1) } else {# otherwise it is just number# so we "just" weighted average them y <-weighted.mean(x, c(237, 187)) y <-round(y, 1) } } else { y <-"\u00b1" }return(y)}
Summary of dataset characteristics for Post Breeding and Post Molting
foraging trips in northern elephant seals.
# seals
All dives
Benthic dives
% Benthic dives per area
Total
w/ loc
w/o loc
Total
Proportion
Shark
Orca
Total
Overlap
Post Breeding
235
1,058,860
1,032,565
26,295
36,443
3.4%
79.8%
100.0%
100.0%
79.8%
Post Molting
187
2,207,361
2,032,528
174,833
44,594
2.0%
90.5%
100.0%
100.0%
90.5%
OVERALL
422
3,266,221
3,065,093
201,128
81,037
2.5%
87%
100%
100%
87%
Code
# display second tabletable_1_b
Table 2:
Individual-level summary of dataset characteristics for Post Breeding
and Post Molting foraging trips in northern elephant seals
# dives
# benthic dives
Hour departure
Hour arrival
Mean
±
SD
Mean
±
SD
Mean
±
SD
Mean
±
SD
Post Breeding
4505.8
±
905.7
155.1
±
275.7
19.0
±
1.8
15.8
±
2.5
Post Molting
11804.1
±
1379.4
238.5
±
471.3
19.2
±
1.7
13.8
±
2.0
OVERALL
7724.6
±
1114.6
191.9
±
362
19.1
±
1.8
14.9
±
2.3
Notes
Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.
Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals
Code
# get the hours of departure and arrivalhour_data <- data_dive %>%group_by(trip, id) %>%summarise(hour_day_departure =first(hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)),hour_day_arrival =last(hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)),.groups ="drop" ) %>%filter(!is.na(hour_day_departure) &!is.na(hour_day_arrival))# compare the circular mean hour of departure between two types of tripwatson.williams.test(list(pb = hour_data %>%filter(trip =="Post Breeding") %>%pull(hour_day_departure) %>%circular(., units ="hours"),pm = hour_data %>%filter(trip =="Post Molting") %>%pull(hour_day_departure) %>%circular(., units ="hours") ))
Watson-Williams test for homogeneity of means
data: pb and pm
F = 0.12556, df1 = 1, df2 = 399, p-value = 0.7233
sample estimates:
Circular Data:
Type = angles
Units = hours
Template = none
Modulo = asis
Zero = 0
Rotation = counter
mean of pb mean of pm
-5.033036 -4.759905
Code
watson.wheeler.test(list(pb = hour_data %>%filter(trip =="Post Breeding") %>%pull(hour_day_departure) %>%circular(., units ="hours"),pm = hour_data %>%filter(trip =="Post Molting") %>%pull(hour_day_departure) %>%circular(., units ="hours") ))
Watson-Wheeler test for homogeneity of angles
data: pb and pm
W = 0.14773, df = 2, p-value = 0.9288
Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.
Code
# to choose an individual ~180 dives for 36 hoursdata_dive %>%# rename date as datetimemutate(datetime = date) %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(datetime), datetime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(datetime), datetime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") )) %>%group_by(id, trip) %>%summarise(nb_dives =length(DiveNumber)) %>%arrange(nb_dives) %>%View()
Code
# import the csv exporttdr_data <-lapply(list.files(path ="../data/",pattern ="*_iknos_raw_data.csv",full.names = T ),function(x) {# extract id id_name <-last(strsplit(strsplit(x, split ="_")[[1]][1],split ="/" )[[1]])# load file dt <- readr::read_csv(x)# add id dt %>%mutate(id =as.numeric(id_name)) }) %>%bind_rows()# choose the lagnb_days <-1.5# compute the datadata_plot <- tdr_data %>%# convert time into date, a Posixmutate(datetime =matlab_to_posix(time),date =as.Date(datetime) ) %>%# add beginning and end date of the tripmerge( ., data_dive %>%group_by(id, trip) %>%summarise(date_low =min(date),date_high =max(date),.groups ="drop" ),by ="id",all.x =TRUE ) %>%# identify what to keepmutate(to_keep =if_else(between(datetime, date_low, date_high), T, F)) %>%# keep only between the beginning and the endfilter(to_keep) %>%# sort to make sure the difftime are being done correctlyarrange(id, time) %>%# merge to get lat and lonmerge( ., data_dive %>%group_by(id, date =as.Date(date)) %>%summarise(lat =median(Lat, na.rm = T),lon =median(Long, na.rm = T),.groups ="drop" ),by =c("id", "date"),all.x =TRUE ) %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(datetime), datetime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(datetime), datetime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") ))# add bathysetDT(data_plot)setDT(data_dive)data_dive[, datetime := date]setkeyv(data_plot, c("id", "datetime"))setkeyv(data_dive, c("id", "datetime"))data_plot <-copy(data_dive) %>% .[, .(id, datetime, bathy, DiveTypeName)] %>% .[data_plot, roll = T]# get sunrise sunset informationsun_data <- data_plot %>%select(date, lat, lon, id, phase) %>%unique() %>%ungroup() %>%mutate(getSunlightTimes(data = ., keep =c("sunrise", "sunset"))) %>%# fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1mutate(sunset =if_else( # if sunset is the next dayabs(as.numeric(difftime(as.POSIXct(date), sunset, units ="days") )) >1,# subtract 1 sunset - (1*60*60*24),# otherwise leave it that way sunset ))# trim if night beging before the trip or end after they arrivesun_data <-merge( sun_data %>%select(date, lat, lon, id, phase, sunrise, sunset), data_plot %>%group_by(id, phase, trip) %>%summarise(date_low =min(datetime),date_high =max(datetime),.groups ="drop" ),by =c("id", "phase"),all.x =TRUE, ) %>%mutate(# make sure sunset is not outside the 36h periodcorrect_sunset =if_else(# if notbetween(sunset, date_low, date_high), sunset,# if so, replace by one extremumif_else(sunset < date_low, date_low, date_high) ),# make sure sunrise is not outside the 36h periodcorrect_sunrise =if_else(# if notbetween(sunrise, date_low, date_high), sunrise,# if so, replace by one extremumif_else(sunrise < date_low, date_low, date_high) ) ) %>%# shift everything to start 2000-01-01 00:00:00mutate(correct_sunset_shift = correct_sunset - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ),correct_sunrise_shift = correct_sunrise - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ) )# plotdeparture_dp <- data_plot %>%filter(phase =="Departure") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(datetime = datetime - (first(datetime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = datetime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Departure"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = datetime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("+6h", "+12h", "+18h", "+24h", "+30h"),position ="top" ) +scale_y_reverse(limits =rev(range(data_plot$depth))) +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free_x",# independent = "x" ) +theme(strip.placement ="outside",strip.text.y =element_blank(),legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",# arrow = arrow(length = unit(0.2, "lines"), type = "closed") ),axis.line.y =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"),type ="closed",ends ="first" ) ),axis.title.x =element_blank() )arrival_dp <- data_plot %>%filter(phase =="Arrival") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(datetime = datetime - (first(datetime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = datetime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Arrival"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = datetime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("-30h", "-24h", "-18h", "-12h", "-6h"),position ="top" ) +scale_y_reverse(limits =rev(range(data_plot$depth))) +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free_x",# independent = "x" ) +theme(strip.placement ="outside",legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"), type ="closed") ),axis.line.y =element_blank(),axis.ticks.y =element_blank(),axis.text.y =element_blank(),axis.title.y =element_blank(),axis.title.x =element_blank() )fig_label_x <-ggplot(data.frame(l ="Time", x =1, y =1)) +geom_text(aes(x, y, label = l)) +theme_void() +coord_cartesian(clip ="off")
Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.
Notes
TODO: Add three dots between Departure and Arrival
4.3 Figure 2
Let’s create a table specific for this figure that only contains for each individual:
the first dive
the last dive
all benthic dives
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# then by individualgroup_by(id) %>%# keep only the first, last date, but also all benthic divesfilter(date ==min(date) | date ==max(date)) %>%# ungroupungroup() %>%# change colnames to make getSunlightTimes worksmutate(date =as.Date(date_tz),lat = Lat,lon = Long ) %>%# calculate hour of sunset and hour of sunrisemutate(getSunlightTimes(data = .,keep =c("sunrise", "sunset"),tz =tz_lookup_coords(lat = Lat,lon = Long,method ="accurate" ) )) %>%# extract hours of sunset an sunrisemutate(h_sunrise =round(hour(sunrise) +minute(sunrise) /60+second(sunrise) / (60*60)),h_sunset =round(hour(sunset) +minute(sunset) /60+second(sunset) / (60*60)) )
Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.
Note
n = 401, because we removed individuals without location data (filter(!is.na(Lat)))
# functio to convert time in radian (probably already exist in some package...)hour_time_2_radian <-function(x, high =24, low =0) {# https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles ptp <- high - low rad <- x * (2* pi) / ptp rad <- (rad + pi) %% (2* pi) - pireturn(rad)}
(Not meant to be in the paper) Average hour of departure and arrival
for post breeding and post molting trip with the result of the
associated Rayleigh test.
Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.
n = 401, because we removed individuals without location data and those that do not have any benthic dives (filter(!is.na(Lat) & DiveTypeName == "Benthic"))
4.5 Figure 4
Code
# https://mycolor.space/?hex=%23008C49&sub=1# build palette colorcols_ind_1 <-c("#008c49", "#00c9d3", "#bffcfd", "#458081")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%2300C9D3-%23BFFCFD-%23458081cols_ind_2 <-c("#008c49", "#002e00", "#00c9d3", "#deb887")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%23002E00-%2300C9D3-%23DEB887
Code
# main plotbathy_prop <- data_dive %>%# keep only rows with location datafilter(!is.na(Lat)) %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# positive bathymutate(bathy =abs(bathy)) %>%# create class of bathymetrymutate(bathy_class =cut( bathy,seq(0, 6000, 400),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per bathy_classgroup_by(bathy_class) %>%# to get the Proportion of different dive types per bathy_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = bathy_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Bathymetry (m)",y ="Dive type proportion (%)" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="top",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # second plot# bathy_hist <- data_dive %>%# # keep only negative bathy# filter(bathy < 0) %>%# # and remove outliers# filter(bathy > -6000) %>%# # create class of bathymetry# mutate(bathy_class = fct_rev(cut(# bathy,# seq(-6000, 0, 400),# ordered_result = T,# dig.lab = 4# ))) %>%# # calculate by bath_class and animal# group_by(bathy_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = bathy_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))
Code
# main plotdist_coast_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 1900, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from the coast (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_coast_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_coast * 1000 * 111 <= 1900) %>%# # create class of bathymetry# mutate(dist_class = cut(# # convert decimal degree/1000# dist_coast * 1000 * 111,# seq(0, 2000, 100),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## the actual plot# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))
Code
# the main plotdist_dep_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_dep >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from departure (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_dep_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_dep / 1000 < 4200) %>%# # create class of bathymetry# mutate(dist_class = cut(# dist_dep / 1000,# seq(0, 4200, 300),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))
# display(guide_area() + bathy_prop + dist_coast_prop + dist_dep_prop) +theme(plot.margin =unit(c(0, .02, -1, .05), units ="npc")) +plot_annotation(tag_levels ="A") +plot_layout(nrow =5, guides ="collect")grid::grid.draw(grid::textGrob(ylab, x =0.07, y =0.55, rot =90))
Figure 4: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).
Notes
TODO: change the color code to match the color we attributed to benthic dives (greens?)
n = 401, because we removed individuals without location data (filter(!is.na(Lat)))
4.6 Figure 5
Code
# the last plotprop_at_sea_dt <- data_dive %>%group_by(id) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure) %>%# perform our calculation per sealgroup_by(id) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class )# plotprop_at_sea <- prop_at_sea_dt %>%# intiate plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +annotate("text", x =5.3, y =0.4, label ="Departure") +annotate("text", x =15, y =0.5, label ="Arrival") +# geom_curve(# data = tibble(# x1 = c(4, 17),# x2 = c(1.5, 19.5),# y1 = c(0.4, 0.5),# y2 = c(0.31, 0.55)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +geom_curve(data =tibble(x1 =c(4),x2 =c(1),y1 =c(0.4),y2 =c(0.33) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =0.3 ) +geom_curve(data =tibble(x1 =c(16),x2 =c(19.7),y1 =c(0.5),y2 =c(0.57) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =-0.56 ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea
Figure 5: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.
Notes
?
n = 422
Test to compare if the proportion of benthic dives at arrival is higher than the one at departure:
\(H_{0}: p_{departure} \ge p_{arrival}\)
\(H_{1}: p_{departure} \lt p_{arrival}\)
Code
# test to compare the proportion between arrival and departurez_test <-prop.test(x =c( prop_at_sea_dt %>%filter(day_class ==first(day_class)) %>%pull(nb_benthic_class), prop_at_sea_dt %>%filter(day_class ==last(day_class)) %>%pull(nb_benthic_class) ),n =c( prop_at_sea_dt %>%filter(day_class ==first(day_class)) %>%pull(nb_total_class), prop_at_sea_dt %>%filter(day_class ==last(day_class)) %>%pull(nb_total_class) ),alternative ="less")z_test
2-sample test for equality of proportions with continuity correction
data: c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_benthic_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_benthic_class)) out of c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_total_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_total_class))
X-squared = 5973, df = 1, p-value < 2.2e-16
alternative hypothesis: less
95 percent confidence interval:
-1.000000 -0.236269
sample estimates:
prop 1 prop 2
0.3093339 0.5505404
Code
# extract the z statisticsqrt(z_test$statistic)
X-squared
77.28513
\(H_{0}\) is rejected, so we can say that the porportion of benthic dives at arrival is significantly superior than departure!
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat)) %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName) %>%# create col to nicely display dive typemutate(origin =paste(DiveTypeName, "dives"))
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# make it's group_by origindf_kernel_dens_sf <-group_by(df_kernel_dens_sf, origin)# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plotto_plot <- trip_zoom_out +# rename axislabs(x ="Longitude", y ="Latitude") +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# facet by originfacet_wrap(. ~ origin)# reorder layersto_plot$layers <-if (with_bathy) { to_plot$layers[c(1, 2, 4, 3)]} else { to_plot$layers[c(1, 3, 2)]}# displayto_plot
Figure 6: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probs represents different levels of probability density.
Notes
lmk!
4.7.3 Figure 2
Code
# the last plotprop_at_sea_trip <- data_dive %>%group_by(id, trip) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id, trip) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure, trip) %>%# perform our calculation per sealgroup_by(id, trip) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class, trip) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +facet_wrap(. ~ trip) +# annotate("text", x = 5.3, y = 0.4, label = "Departure") +# annotate("text", x = 15, y = 0.5, label = "Arrival") +# geom_curve(# data = tibble(# x1 = c(4),# x2 = c(1),# y1 = c(0.4),# y2 = c(0.33)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +# geom_curve(# data = tibble(# x1 = c(16),# x2 = c(19.7),# y1 = c(0.5),# y2 = c(0.57)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = -0.56# ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea_trip
Figure 7: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea during post breeding and post molting trip. Time spent at sea is classified into intervals of five percent for graphical representation.